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Pairing Correlations in the two-layer attractive Hubbard Model 
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Studies of systems with two fermionic bands with repulsive interaction strength U have a long history, with 
the Periodic Anderson Model (PAM) being one of the most frequently considered Hamiltonians. In this paper, 
we use Quantum Monte Carlo to study analogous issues for attractive interactions. As in the Periodic Anderson 
Model, we focus on a case where one band is uncorrelated (U = 0), and focus on the effect of hybridization V 
between the bands on the pairing correlations. A key difference with the PAM is that there is no sign problem, 
so that we are able to explore the physics of doped multi-band attractive systems at low temperatures whereas 
ground state properties of repulsive models can be determined only at half-filling. For small V, pairing in the 
U < layer induces pairing in the [7 = layer. At larger V the ground state of the coupled system loses 
its superconducting character. The Quantum Monte Carlo data are complemented by results obtained with the 
Bogoliubov-de Gennes approximation. 

PACS numbers: 74.20.-z,71.10.Pm,74.25.Dw,74.78.-w 



I. INTRODUCTION 

Cuprate superconductors are characterized by CUO2 planes 
which are rather isolated from each other by intervening rare 
earth atoms. The larger ("c-axis") separation of Cu atoms 
perpendicular to the planes compared to the in-plane lat- 
tice constants (a, b) has focused theoretical attention on the 
magnetic and superconducting properties of two-dimensional 
models, most notably the square lattice Heisenberg and Hub- 
bard Hamiltonians. Indeed, one of the most fundamental the- 
oretical questions which arose in the initial investigations of 
possible models of high temperature superconductivity con- 
cerned the nature of the magnetism in the ground states of 
these Hamiltonians: Were they spin-liquids or did they ex- 
hibit long range order? It was established numerically™ 
that both the Heisenberg Hamiltonian and the half-filled Hub- 
bard Hamiltonian have long range antiferromagnetic order at 
T = 0. A still-unresolved question is the nature of pairing in 
the doped Hubbard Hamiltonian. 

A natural extension of these single layer questions con- 
cerns the behavior of spin and pairing correlations in coupled 
planes. Interlayer connections are required to elevate the mag- 
netic ordering, which is possible only at T = in two di- 
mensions, to the finite temperature Neel transitions observed 
experimentally. Similarly, the superconducting transition tem- 
perature tends to increase with the number of adjacent CUO2 
layers between the charge resevoirs, getting steadily larger 
from single layered La2_ K Sr x Cu04, to bilayer YBaCu^Oe+y 
to the HgBaCaCuO sequence where T c peaks at about 130° 
K for materials with n = 3 layers^. Antiferromagnetism 
and thus superconductivity exhibit the expected behavior that 



adding additional (non-frustrating) couplings, and increasing 
the dimensionality, enhances the tendency to order. 

The route from 2D to 3D is, however, not completely 
straightforward. Quantum Monte Carlo (QMC) studies of 
coupled Heisenberg layers*, have shown that when the inter- 
plane coupling grows sufficiently large, there is a Quantum 
Phase Transition (QPT) from a ground state with long range 
magnetic order to a spin-gapped ground state, in which sin- 
glets form between pairs of spins in the two planes4 A sim- 
ilar effect is seen in QMC studies of the two-layer half-filled 
Hubbard Hamiltonian^. 

The basic phenomenon of singlet formation usurping mag- 
netic order in itinerant fermion Hamiltonians does not cru- 
cially depend on having the same interaction strength U, or 
hopping t, for the two species. Indeed, a Hamiltonian with 
U = in one (conduction) band and hopping t — in an- 
other (localized) band, is the commonly encountered Periodic 
Anderson Model (PAM). At half-filling, for weak hybridiza- 
tion between the conduction and localized bands, antiferro- 
magnetic order is present. As the hybridization grows big- 
ger, however, "Kondo" singlets form and destroy antiferro- 
magnetic (AF) order. Thus the PAM exhibits a similar QPT 
as for two identical Hubbard layers (t\ — t-x and U\ = Uz). 
This close similarity of the PAM and Hubbard bilayers em- 
phasizes that the physics of fermions in multiple layers can 
equivalently be interpreted in in terms of multiple orbital s. 

The purpose of this paper is to study the possibility of anal- 
ogous phenomena in the case when there is an attractive in- 
teraction. More specifically, we examine the nature of pairing 
correlations for coupled planes (orbitals) in situations when 
an attraction is present only in one plane (orbital). We use the 
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—U Hubbard Hamiltonian as an appropriate simple model, 
and employ a combination of Quantum Monte Carlo (QMC) 
and Bogoliubov-de Gennes (BdG) Mean Field Theory. We 
will focus on a case analogous to that of the PAM, namely 
when U = in one of the layers. We are interested in both 
how the interlayer hopping affects the s-wave pairing corre- 
lations in a layer with U < 0, and also whether pairing can 
be induced by a "proximity effect" in the U = layer. Such 
behavior would be directly analogous to the "RKKY" spin 
polarization cloud which is induced in the conduction band 
by the presence of the local moments in the PAM. A partic- 
ularly interesting issue is whether the QPT which occurs at 
commensurate filling also occurs in the doped case. 

The remainder of this paper is organized as follows: In 
Sec. II we write down the precise model and give an overview 
of the QMC and BdG computational methodologies. Sees. Ill 
and IV present the BdG and QMC results respectively. Sec. V 
recaps our conclusions. 



H. MODEL AND COMPUTATIONAL METHODS 

Our starting point is the two layer (orbital) attractive Hub- 
bard Hamiltonian, 

<ij>la 

V/^( c ila C i2a + c i2(T C ilcr) — E ^l C ila C ila 
ia ila 

- ENKit-^Ku-^) • (i) 

i 

Here an intralayer kinetic energy term c| Jo .c.j ff describes the 
creation of a fermion with spin a on site i of layer I and its 
destruction on site j of the same layer. A second interlayer 
kinetic energy term c\ la c iVa hybridizes the two layers 1,1'. 
Fermions of different spin feel an attractive interaction — \Ui\ 
on layer I. The chemical potentials fxi control the filling. Our 
lattice geometry consists of two coupled square lattice of lin- 
ear size L. We choose U\ < in the 'superconducting' plane 
and U2 = in the metallic one. 

At half-filling, a particle-hole transformation formalizes the 
similarity of the repulsive and attractive models. The vanish- 
ing of antiferromagnetic order with increasing V for U > 
bilayers immediately implies that ground state pairing (and 
charge density wave order, with which it is degenerate) must 
be destroyed as V grows for U < 0. However, whereas 
AF order occurs only at half-filling (and only at T = 0) 
in the 2D repulsive Hubbard Hamiltonian, superconductiv- 
ity appears away from half-filling, and at a finite (Kosterlitz- 
Thouless) transition temperatur e - 24 ' 25 , in the attractive Hub- 
bard Hamiltonian. This opens up fundamentally new issues 
in the attractive case. Whether the fermions on sites in adja- 
cent layers will still lock into a local object and destroy long 
range order when the filling is incommensurate is an interest- 
ing question. The absence of a sign problem in the case of 
the attractive Hubbard model enables the study of low tem- 



peratures even when the filling is incommensurate, allowing 
access to ground state properties. 

We perform "determinant" QMC 9 simulations of the 
Hamiltonian, Eq.Q]by writing down a path integral for the par- 
tition function Z. The exponential of the interaction term Ui 
is decoupled with a Hubbard-Stratonovich field which allows 
the trace of the remaining exponentials of quadratic forms of 
the fermion operators to be performed analytically. The re- 
sult is a sum over configurations of the Hubbard-Stratonovich 
field with a weight which takes the form of the product of the 
determinants of two matrices of dimension the spatial lattice 
size, one for each spin species. In the case of attractive in- 
teraction the two matrices, and hence their determinants, are 
identical, so that the weight is a perfect square and there is no 
sign problem. 

We present results here for the real space pair correlation 
functions in the two orbitals, 

PiU) = (&i +j ,Ah) 



and also their associated structure factors, 

3 

as a functions of interband hybridization V, on-site attraction 
Ui, density p, and temperature T. The long distance behavior 
of P ; (j) yields the square of the superconducting order pa- 
rameter, as does the lattice size dependence of the structure 
factor. 

Eq.Q]can also be studied in Mean Field Theory via the solu- 
tion of the BdG equations. In this approach, the four fermion 
(interaction) term is decoupled in the Hamiltonian itself, 

U eB = ~ t (4la C jla + C )la C ila) 

<ij>l<7 

V Z^( C ila C i2a + C i2a C ila) ~ E ^il C lla C ila 
ia ila 

- EN[ A *i4t c lu + A* lC?uCllt ] . (4) 

i 

Here the gap A tl = (c ilt c ui ) and density (n u<7 ) = 

(c ila .cua) are determined self-consistently by diagonalizing 
the quadratic BdG Hamiltonian and putting together the 
eigenvalues and eigenvectors appropriately to compute refined 
values which are inserted back in the Hamiltonian in an iter- 
ative process. A related BdG treatment of superconductivity 
in the presence of randomness (spatially varying \xi) is given 
in lllOill ill . p,u = (j, + \Ui\(nu)/2 includes a site-dependent 
Hartree shift with (n«) = ^2 a {nua)- 

In the general case, when inhomogeneous terms are present 
in the Hamiltonian, or are expected to develop spontaneously, 
(like the charge and spin stripes of the doped, repulsive Hub- 
bard model), the order parameter and densities are allowed to 
depend on site index and the diagonalization must be done 
numerically. If translation invariance is present, Eq. |4]can be 
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diagonalized analytically by going to momentum space. The 
resulting momentum sums are typically still done numerically, 
but larger lattices can be studied. We have implemented both 
approaches in the work presented here. 

Although inhomogeneous Hartree-Fock theory allows for 
spatial variation of the densities and superconducting order 
parameters, the expectation values in Eq.|4]are independent of 
imaginary time, unlike the fluctuating Hubbard-Stratonovich 
field in the QMC approach. This leads to a less accurate 
treatment of interparticle correlations, but a numerically much 
more simple problem. In particular, larger spatial lattices can 
be studied, and the BdG solution can exhibit broken symmetry 
so that A itself can be nonzero, rather than having to be ex- 
tracted from the asymptotics of the correlation function Eq. 2. 
Here the combination of QMC and BdG will serve to provide 
complementary information. 

In the case of the BdG calculations our focus will be on 
the value of the order parameter A as a function of the same 
parameters as varied in the QMC. As mentioned above, the 
asymptotic value of Pi(j) is a measure of A 2 . The BdG ap- 
proach, since it neglects fluctuations, yields larger values of A 
than those of the QMC, which is an exact method. 



m. BOGOLIUBOV-DE GENNES RESULTS 

The results of BdG calculations are summarized in Fig. Q] 
We consider two layers with U\/t = —8 and U2 = 0, at a 
density of p = 0.8 electrons per site and study the evolution 
of the pairing correlation, A;, in the two layers as a func- 
tion inter-layer coupling, V. Pairing in the correlated layer is 
suppressed by coupling to the non-correlated layer and even- 
tually destroyed completely at a critical V c , analogous to the 
destruction of antiferromagnetism by singlet formation in the 
repulsive case. A further interesting feature of Fig. [T]is that the 
I = 2 layer, with U2 = 0, develops induced pairing with the 
onset of inter-layer coupling. The pairing amplitude, P2, in- 
creases with increasing V up to a maximum before decreasing 
and eventually vanishing at V c . For V > V c , the ground state 
consists of weakly interacting inter-layer dimers. Crucially, 
the pairing order parameter in both layers remains non-zero 
over a finite range of inter-layer coupling - the induced pairing 
reported previously at half-filling extends to finite doping s. In 
the next section, we present extensive QMC results to confirm 
and complement the BdG results. 



IV. QUANTUM MONTE CARLO RESULTS 

A. Half Filling 

We begin our discussion of QMC results at half-filling, 
px = p2 = 1. As discussed previously, this corresponds to 
the particle-hole symmetric point /.ti = /i2 = where there 
exists an exact mapping to the repulsive Hubbard model and 
the results can be benchmarked against previous studies of the 
Hubbard model [7 8]. Lattices of the form 2 x L x L, with 
4 < L < 14 were studied with periodic boundary conditions 




FIG. 1: Bogoliubov-de Gennes result for 8 x 8 and 12 x 12 bilayer 
lattices, j3 = 8, Ui = —8, p = 0.8. The figure shows that there is 
an induced SC in the noninteracting layer, with negligible finite size 
effects. 



over a wide range of inter-layer hybridization, and on-site 
interaction for the interacting layer. An inverse temperature 
fit = L was found to be sufficient for the observables to have 
converged to their ground state values. 

The spatial dependence of the pair correlations Pi(j) at 
half-filling (p = 1), intersheet hybridization V = 0.75t, and 
interactions U\ = — 6,^/2 = Ois shown in Fig. [2] The sep- 
aration j follows a trajectory along the x axis to maximal a; 
separation (j,0) on a lattice with periodic boundary condi- 
tions, and then to (4, before returning to separation (0, 0). 
Results for lattices with L = 8, 12 are shown in the figure. 
The interacting layer I = 1 exhibits clear long range order 
with Pi (j) nearly independent of j beyond a separation of ap- 
proximately j . The noninteracting layer I = 2 also exhibits 
proximity-effect induced long range order, although the corre- 
lation function is roughly an order of magnitude smaller. This 
corresponds to an order parameter Aj which is approximately 
a factor of three smaller. 

Increasing the inter-layer hybridization reveals that the pair 
structure factor for the interacting layer, , decreases mono- 
tonically and eventually the intra-layer pairing is destroyed at 
some critical V c . On the other hand, the pairing structure fac- 
tor for the non-interacting layer P s 2 varies non-monotonically 
with V. It is vanishingly small at small V. As V increases, 
long-range pairing correlations increase, reach a maximum at 
an intermediate V (which depends on the strength of on-site 
interaction in the interacting layer) and then decrease continu- 
ously to zero at V c . For V > V c , the ground state is dominated 
by inter-layer singlets. Fig.[3]shows the finite size dependence 
of the pairing structure factor for the two layers at several val- 
ues of V. For V > 1.6, both P} and P 2 extrapolate to zero in 
the thermodynamic limit indicating the absence of intra-layer 
pairing in this regime. 

Here and further, the finite size scaling was done by linear 
fit in the plane 1/L — Pj. The scatter of the data in Fig.[3]for 
different values of the interlayer hybridization V is dominated 
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(0,0) (L/2,0) (L/2.L/2) (0,0) 

j=(dx,dy) 

FIG. 2: The dependence of the ground state pair correlation func- 
tions, Pi,2(j), on separation j in a bilayer system comprised of one 
correlated layer coupled to an uncorrelated layer for two different lat- 
tice sizes L. The system is at half-filling, p = 1, with interplane hy- 
bridization V — 0.8t, and on-site interaction strength Ui = —6t for 
the correlated layer. The correlation functions converge to non-zero 
value at large separations, providing clear evidence for long-range 
order. Despite the absence of interactions, there is proximity-effect 
induced long-range order in the uncorrelated layer, / = 2, with an 
order parameter A2 approximately a factor of three lower than for 
the correlated layer. As in Fig. [I] finite size effects are modest. 



by the error bars associated with the linear fit to the structure 
factor on different lattice sizes used in the extrapolation, as 
opposed to the statistical errors for an individual run. The 
data scatter serves as a stand-in to provide a visual indication 
of the uncertainty in the the calculation. 




FIG. 3: The pair structure factor for a bilayer at half-filling for attrac- 
tive on-site interaction in the correlated layer U — — Wt. The main 
panels show the extrapolated structure factor as a function of the 
inter-layer coupling, V. The insets show the extrapolation of the data 
from finite-sized systems to the thermodynamic limit. Top panel: 
The pair structure factor in the correlated layer, P 3 . For inter-layer 
hybridization V < 1.6 the extrapolated structure factor is non-zero 
indicating the existence of long-range order (LRO). Bottom panel: 
The pairing structure factor P s 2 in the non-interacting layer. P 3 in- 
creases from zero for V > 0, signaling the onset of induced pairing. 
With further increase of V, LRO vanishes simultaneously in both 
layers. 



B. Incommensurate Filling 

We now turn to the case of doped planes. As mentioned 
in the Introduction, this is interesting for two reasons. First, 
there is a possibility of a finite temperature KT transition be- 
cause the CDW-SC symmetry is broken. Second, simulations 
of coupled, doped repulsive layers at low T are not possible. 
Hence there is no DQMC information available for the nature 
of magnetism in coupled layers away from half-filling. Our 
simulations address this issue, albeit for attractive on-site in- 
teractions. 

The results in Fig. [4] confirm unambiguously that induced 
pairing in the uncorrelated layer extends to finite dopings 
away from half-filling; indeed the pair correlation functions 
in the two layers are qualitatively similar to those at half- 
filling. Pairing in the correlated layer decreases monotoni- 
cally to zero at a finite V c and that in the non-interacting layer 
increases rapidly from zero as the inter-layer hybridization is 
turned on, reaches a maximum value at intermediate V and 
then decreases to zero. At intermediate to strong values of 
the on-site interaction (\U\ > 6t), the ground state of a sys- 
tem of coupled correlated and uncorrelated layers away from 



half-filling consists of intra-layer pair formation in both layers 
over a finite non-zero range of inter-layer coupling. Eventu- 
ally, at sufficiently strong inter-layer hybridization, pairing in 
both layers is destroyed and the ground state is dominated by 
singlet formation between the layers. Induced pairing is found 
to be absent for \U\ < 6t. The simulation results are summa- 
rized in Fig. [4] where the extrapolated pairing structure factor 
for the two layers are shown as a function of the inter-layer 
hybridization, V, for some representative values of the on-site 
interaction in the correlated layer at a fixed value of the den- 
sity, p = 0.8. The maximum of the pairing structure factor 
in the uncorrelated layer and the critical value of inter-layer 
coupling for complete suppression of pairing, V c , increase 
with \U\ . Our results are consistent with the disappearance 
of pairing in both the layers occurring simultaneously at V c . 
The error associated with determining the extrapolated pair- 
ing structure factor in the thermodynamic limit limits the ac- 
curacy with which V c can be determined to about ±0.2i. In a 
marked departure from the BdG results, the QMC data shows 
unambiguously that a finite non-zero inter-layer coupling is 
required for the onset of induced pairing in the uncorrelated 
layer. This is a consequence of the mean field nature of the 
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FIG. 4: The pair structure factor for a bilayer away from half -filling, 
p = 0.8, for three representative values of the attractive on-site inter- 
action in the correlated layer: Ui = —At, —6t, — lOt. For \Ui | < 6t, 
Ps = at all values of V, indicating the absence of induced pairing. 
At stronger on-site interaction (|t/i| > 6t) in the correlated layer, 
induced pairing appears in the uncorrelated layer over a finite range 
of V. The critical V c for destruction of intra-layer pairing increases 
with \Ui\. Insets: Finite size scaling for Ui = — lOt. 



BdG results. 

Fig.|4]shows that Pi increases systematically as U changes 
from U = —4 to U = — 10 and, indeed, shows no sign of sat- 
uration. This is at variance with the behavior of the magnetic 
response of the repulsive Hubbard Hamiltonian which first in- 
creases with U but then reaches a maximum at U ~ 8t before 
falling. This behavior is understood qualitatively from the fact 
that the superexchange J = 4t 2 /U which provides the large 
U magnetic energy scale declines with U. Thus, for example 
the Neel temperature of the 3D Hubbard model at half-filling 
is maximized at U ~ 8t. The analog in the attractive Hub- 
bard model is the fact that the pairs become heavy, with an 
effective hopping t g ~ t 2 /U associated with the fact that 
a pair must be (temporarily) broken before it can hop. This 
analogy suggests the pairing response might also be maximal 
at an intermediate U. However, because of the sign problem, 
it is not known from DQMC whether such a non-monotonic 
behavior exists in the repulsive model away from half-filling. 
Our results here for the attractive case, which show no sign of 
saturation with increasing U, suggest that it might not. 

Having established the existence of induced pair formation 
at finite doping and the evolution of the associated structure 
factor with V, we investigate the dependence of the pairing 
phenomenon on the density of electrons in the layers. Our re- 
sults (Fig. [5]) show that induced pairing occurs generically at 
all densities. For all density values studied, the pairing strac- 



FIG. 5: The pair structure factors at various filling values. Induced 
pairing occurs at all densities, with the amplitude of induced pairing 
increasing with doping away from half-filling. Insets: Finite size 
scaling for p = 0.8. 



ture factor behaves in a qualitatively similar manner, with the 
pairing strength optimally enhanced at p = 0.8, as evidenced 
by the initial increase in peak height with decreasing density. 
Pairing appears to show a modest decline as the density is re- 
duced further to p — 0.6, although the effect is rather small 
given the scatter in the data. Such a non-monotonic doping 
dependence would be consistent with the behavior of the su- 
perconducting T c with doping for a single 2D sheet. T c rises 
abruptly from zero at half-filling where charge density order 
competes with pairing, but then shows a gradual decline be- 
low an optimal doping^. 



V. CONCLUSIONS 

In this paper we have used Quantum Monte Carlo and 
Bogoliubov-de Gennes mean field theory to elucidate the 
properties of a multi-orbital attractive Hubbard Hamiltonian. 
Our key conclusions are the quantification of the induced pair- 
ing in a noninteracting orbital by superconductivity in a cor- 
related orbital, and, conversely, the suppression of pairing in 
an orbital with an attractive interaction through its coupling to 
a noninteracting orbital. 

As we have emphasized in the introduction, analogous is- 
sues for repulsive models have a long history. In fact, over 
the last 4-5 years, the question of 'orbitally selective Mott 
transitions', has been extensively explored. The fundamental 
objective has been an examination of the density of states to 
determine whether or not one subsystem can be metallic (zero 
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energy gap at E-p) and the other insulating (nonzero gap), or 
whether the coupling forces the transition to occur simulta- 
neously in all orbitals (layers)^— . Our work suggested that 
for attractive interactions the destruction of superconductivity 
at large inter-orbital hybridization V always occurs simultane- 
ously. It is to be noted, however, that it appears that at weak V 
pairing can exist in the correlated orbital before it is induced 
in the noninteracting one. That is, we find that for the attrac- 
tive case at small interlayer coupling, SC exists only in the 
interacting layer despite the fact that at larger V the vanishing 
of SC appears at a common V c . 

Our work is related to several other recent determinant 
QMC and Lanczos studies of the attractive Hubbard model. 
Assman et.d@& found that in an attractive Hubbard model 
with a smoothly varying chemical potential (such as would be 
present in a confined ultracold atomic cloud) pairing is signifi- 
cantly increased in the half-filled portion of the lattice through 
its contact with doped regions. That is, the suppression of 
superconductivity caused by the appearance of a degenerate 



charge density wave phase at p = 1 is eliminated. Paiva et.al. 
explored a one-dimensional model of borocarbides in which 
U < and [7 = sites alternate 2 ?. Superconductivity is 
found to be possible only above a critical density. 

Finally we note that it was recently shown in repulsive mod- 
els with more than two layers that there can be further interest- 
ing and even surprising evolution of magnetic properties with 
the hopping across an interface between a Mott insulator and a 
metal 2 ^. It would be interesting to pursue related questions in 
the repulsive model and, in particular whether the formation 
of a layer of local pairs at the boundary can lead to a shielding 
of penetration effects. 

We thank Simone Chiesa and Axel Euverte for very help- 
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Award W911NF07 10576 with funds from the DARPA OLE 
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search grant; and by the CNRS-National University of Sin- 
gapore FSQL LIA joint research grant. 
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